function [x,output]=agentsolverho_4agents(xguess,index,ewsb,waut,bind,nobind);
global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table
objval = [];
output=[];
options = optimset('Display','iter','TolFun',1e-8,'GradObj','on');
[x,fval,exitflag]=fsolve(@objfun,xguess,options);
output(6)=exitflag;
function f = objfun(x);
i=index(1);
j=index(2);
consu(1:agent-1) = x(1:agent-1);
consu(agent)=Y(j)-sum(x(1:agent-1));
f=[];
wint=interpolation_4agents(ewsb, index,x./Y(j));
for g=1:agent
output(g)=utility(rho,consu(g))+beta*wint(g);
end
for k=1:length(nobind)
f = [f; gammagrid(i,nobind(end))/gammagrid(i,nobind(k))-consu(nobind(end))/consu(nobind(k))];
end
for k=1:length(bind)
f = [f; ((utility(rho,S(j,bind(k)))+beta*P(j,:)*waut(:,bind(k))) - output(bind(k)))];
end


end

end